# ===============================================================#
#                     Replication files for:                     #
#.  "Attitudinal and Behavioral Legacies of Wartime Violence:    #
#                      A Meta-Analysis"                          #
#                        Joan Barceló                            #
#               American Political Science Review                #
#               Last update: September 3, 2025                   #
# ===============================================================#

#################################
# Appendix G.2 / Table G.1: Egger's Test for Publication Bias
################################

## ------------------------ Load Packages ------------------------

library(dplyr)
library(tidyr)
library(metafor)
library(clubSandwich)
library(knitr)
library(kableExtra)

# ---------- Settings ----------
min_n <- 10  # minimum number of UNIQUE studies required

# ---------- Paths and file maps ----------
BASE <- "~/Datasets/"

read_meta <- function(path) read.csv(file.path(BASE, path), sep = ",", stringsAsFactors = FALSE)

# Outcome -> filename maps

files_full <- c(
  "Social groups participation" = "meta.groups.csv",
  "Voting"                       = "meta.voting.csv",
  "Generalized trust"            = "meta.generalized_trust.csv",
  "Political participation"      = "meta.political_participation.csv",
  "Interest / knowledge"         = "meta.political_interest.csv",
  "Normative prosociality"       = "meta.normative.csv",
  "Community leadership"         = "meta.leadership.csv",
  "Altruism"                     = "meta.altruism.csv",
  "Attitudes toward wartime enemies" = "meta.warenemies.csv",
  "Threat perceptions"           = "meta.threat.csv",
  "Intergroup attitudes"         = "meta.intergroup.csv",
  "Ingroup trust"                = "meta.ingroup_trust.csv",
  "Group voting"                 = "meta.groupvoting.csv",
  "Group identification"         = "meta.groupid.csv",
  "Political intolerance"        = "meta.political_intolerance.csv",
  "Authoritarian attitudes"      = "meta.authoritarian.csv",
  "Hawkish security preferences" = "meta.hawkish.csv",
  "Support for punitive justice" = "meta.punitive.csv",
  "Extreme ideology"             = "meta.xtrideology.csv",
  "Social intolerance"           = "meta.social_intolerance.csv",
  "Antagonism toward peace"      = "meta.antipeace.csv",
  "Institutional distrust"       = "meta.instmistrust.csv"
)

files_biv <- sub("\\.csv$", "_biv.csv", files_full)
files_qe  <- sub("\\.csv$", "_random.csv", files_full)
files_qe  <- files_qe[!names(files_qe) %in% c("Normative prosociality", "Political intolerance")]

# Load datasets
datasets_full <- lapply(files_full, read_meta)
datasets_biv  <- lapply(files_biv,  read_meta)
datasets_qe   <- lapply(files_qe,   read_meta)

prep_meta <- function(df) {
  if ("coef" %in% names(df) && "se" %in% names(df)) {
    df <- df %>% mutate(yi = as.numeric(coef), vi = as.numeric(se)^2)
  }
  if (!("yi" %in% names(df) && "vi" %in% names(df))) return(NULL)
  if (!"study_id" %in% names(df)) {
    df$study_id <- if ("authoryear" %in% names(df)) df$authoryear else paste0("Study_", seq_len(nrow(df)))
  }
  df$effect_id <- seq_len(nrow(df))
  df
}

egger_one <- function(df, min_n = 10) {
  df <- prep_meta(df)
  if (is.null(df)) return(c(estimate = NA, se = NA, p = NA))
  
  # keep only valid rows for counting and regression
  valid <- is.finite(df$yi) & is.finite(df$vi) & df$vi > 0
  if (!any(valid)) return(c(estimate = NA, se = NA, p = NA))
  df2 <- df[valid, , drop = FALSE]
  
  # count UNIQUE studies using authoryear when available
  if ("authoryear" %in% names(df2)) {
    k_studies <- length(unique(as.character(df2$authoryear)))
  } else {
    k_studies <- length(unique(as.character(df2$study_id)))
  }
  if (k_studies < min_n) return(c(estimate = NA, se = NA, p = NA))
  
  df2 <- df2 %>% mutate(z = yi / sqrt(vi), precision = 1 / sqrt(vi))
  
  lm_fit <- lm(z ~ precision, data = df2)
  rob <- clubSandwich::coef_test(lm_fit, cluster = df2$study_id, vcov = "CR2", test = "Satterthwaite")
  row <- rob[rob$Coef == "(Intercept)", , drop = FALSE]
  
  c(estimate = unname(row$beta), se = unname(row$SE), p = unname(row$p_Satt))
}

fmt_cell <- function(est, se, p) {
  if (is.na(est) || is.na(se) || is.na(p)) return("-")
  stars <- ifelse(p < .001, "***", ifelse(p < .01, "**", ifelse(p < .05, "*", "")))
  sprintf("%.2f%s (%.2f)", est, stars, se)
}

# -------- Run Egger across sets --------
run_set <- function(dset) {
  outs <- names(dset)
  res <- lapply(dset, egger_one, min_n = min_n)
  mat <- do.call(rbind, res)
  as.data.frame(mat) %>% mutate(Outcome = outs, .before = 1)
}

full_res <- run_set(datasets_full)
biv_res  <- run_set(datasets_biv)
qe_res   <- run_set(datasets_qe)

# Merge into one table of formatted cells
tab <- full_res %>%
  select(Outcome, estimate_f = estimate, se_f = se, p_f = p) %>%
  left_join(biv_res %>% select(Outcome, estimate_b = estimate, se_b = se, p_b = p), by = "Outcome") %>%
  left_join(qe_res  %>% select(Outcome, estimate_q = estimate, se_q = se, p_q = p), by = "Outcome") %>%
  mutate(
    Full = mapply(fmt_cell, estimate_f, se_f, p_f),
    Bivariate = mapply(fmt_cell, estimate_b, se_b, p_b),
    `Quasi-Experimental` = mapply(fmt_cell, estimate_q, se_q, p_q)
  ) %>%
  select(Outcome, Full, Bivariate, `Quasi-Experimental`)

# Order rows and add section breaks like the figure
order_vec <- c(
  "Social groups participation","Voting","Generalized trust","Political participation",
  "Interest / knowledge","Normative prosociality","Community leadership","Altruism",
  "Attitudes toward wartime enemies","Threat perceptions","Intergroup attitudes",
  "Ingroup trust","Group voting","Group identification",
  "Political intolerance","Authoritarian attitudes","Hawkish security preferences",
  "Support for punitive justice","Extreme ideology","Social intolerance",
  "Antagonism toward peace","Institutional distrust"
)
tab <- tab %>% slice(match(order_vec, Outcome))

# Fisher omnibus tests per column
fisher_fmt <- function(pvec) {
  pvec <- suppressWarnings(as.numeric(pvec))
  pvec <- pvec[is.finite(pvec) & !is.na(pvec) & pvec > 0]
  if (!length(pvec)) return("-")
  stat <- -2 * sum(log(pvec))
  df <- 2 * length(pvec)
  paste0("\u03C7^2 = ", round(stat, 2), ", df = ", df, ", p = ",
         signif(pchisq(stat, df = df, lower.tail = FALSE), 3))
}

f_full <- fisher_fmt(full_res$p)
f_biv  <- fisher_fmt(biv_res$p)
f_qe   <- fisher_fmt(qe_res$p)

tab_out <- bind_rows(
  tab,
  tibble(Outcome = "Fisher's omnibus test",
         Full = f_full, Bivariate = f_biv, `Quasi-Experimental` = f_qe)
)

print(tab_out)

